#!/bin/bash

set -o errexit
export LC_ALL=C


export BASE_DATA_PATH=/home/mixagol/data/
export BLAST_DB_PATH="${BASE_DATA_PATH}/1_databases/blast_db/GENOMES_DB"

export NORM_MATRIX_DIR=${BASE_DATA_PATH}/3_normed_matrix

export LC_ALL=C
export GENOMES_NUM=`cat "${BASE_DATA_PATH}/3_normed_matrix/map_int_genom.txt" | wc -l`
#g++ -O3 6_annotation/intersect_genes.cpp -o 6_annotation/intersect_genes
#g++ -O3 6_annotation/intersect_genes_avg.cpp -o 6_annotation/intersect_genes_avg

export ANNOTATE_DIR=$1/annot_$(date +%s.%N)
mkdir -p $ANNOTATE_DIR

if [ $2 ]; then
    cat <(echo ">"$2) /dev/stdin > $ANNOTATE_DIR/gene.faa
else
    cat /dev/stdin > $ANNOTATE_DIR/gene.faa
fi

echo -e $(cat $ANNOTATE_DIR/gene.faa | head -n1 | cut -b2-)"\t"$(cat $ANNOTATE_DIR/gene.faa | head -n2 | tail -n1) \
    > $ANNOTATE_DIR/gene.txt

cat $ANNOTATE_DIR/gene.txt | tools/create_bsddb.py $ANNOTATE_DIR/gen.db

cat $ANNOTATE_DIR/gene.faa | 2_raw_matrix_blast/run_blast.sh $BLAST_DB_PATH > $ANNOTATE_DIR/blast_results.txt

#for MNT in {01..51}; do
#    export E_VALUE_RECOMMENDED=1e-"$(echo $MNT | awk '{printf("%02d\n", $1)}')"
#    VSIZE=$(cat $ANNOTATE_DIR/blast_results.txt | awk -F'\t' '{if ($11<='$E_VALUE_RECOMMENDED') print $1"\t"$2}' \
#        | 3_normed_matrix/apply_maps.py \
#             --genomes-map ${BASE_DATA_PATH}/3_normed_matrix/map_int_genom.txt \
#        | 3_normed_matrix/create_matrix.py \
#        | cut -f2 \
#        | awk '{print split($1, M, ",")}')
#
#    export MNT_RECOMMENDED=$MNT
#    if [ $VSIZE -lt 25 ]; then
#        break
#    fi
#done

#for MNT in {{01..15},20,30,40,50}; do
#    if [ $MNT -gt $MNT_RECOMMENDED ]; then
#        continue
#    fi
    export E_VALUE="1e-05"
    export E_DIR=$ANNOTATE_DIR/e_$E_VALUE
    mkdir -p $E_DIR

    cat $ANNOTATE_DIR/blast_results.txt | awk -F'\t' '{if ($11<='$E_VALUE') print $1"\t"$2}' \
        | 3_normed_matrix/apply_genomes_map.py \
            --genomes-map ${NORM_MATRIX_DIR}/map_int_genom.txt \
        | 3_normed_matrix/create_matrix.py \
        > $E_DIR/vectors.txt

    #if [ $(cut -f2 $E_DIR/vectors.txt|awk '{print split($1, M, ",")}') -eq 0 ]; then
    #    rm -r $E_DIR
    #    break
    #fi
    
    6_annotation/intersect_genes \
        -m ${NORM_MATRIX_DIR}/matrix_${E_VALUE}_.txt \
        -v $E_DIR/vectors.txt \
        -g $GENOMES_NUM \
        | gzip \
        > $E_DIR/joined_raw.txt.gz
    
    join -t$'\t' -1 1 -2 2 \
        <(cat ${NORM_MATRIX_DIR}/map_int_gene.txt | sort -t$'\t' -k1,1 -S100M) \
        <(zcat ${E_DIR}/joined_raw.txt.gz | sort -t$'\t' -k2,2 -S100M) \
        | cut -f2,4 \
        | sort -t$'\t' -k1,1 -S100M \
        | gzip \
        > ${E_DIR}/joined.txt.gz
    
    export PROB_LIMIT=1e-7
    
    join -t$'\t' \
        <(zcat ${E_DIR}/joined.txt.gz) \
        ${BASE_DATA_PATH}/5_pathways/gene_kegg_map.txt \
        | awk -F'\t' '{if ($2<'"$PROB_LIMIT"') print $0}' \
        | gzip \
        > ${E_DIR}/joined_kegg_limited.txt.gz
    
    ##
    ## probability<TAB><gene_kegg_id<TAB>function
    ##
    join -t$'\t' \
        <(zcat ${E_DIR}/joined_kegg_limited.txt.gz) \
        <(cat ${BASE_DATA_PATH}/1_databases/genes_info.txt ) \
        | cut -f2,3,6 \
        | sort -t$'\t' -g -k1,1 -S100M \
        | gzip \
        > ${E_DIR}/results_gene_functions.txt.gz
   
    ##
    ## probability<TAB><gene_kegg_id<TAB>function
    ##
    join -t$'\t' \
        <(zcat ${E_DIR}/joined_kegg_limited.txt.gz) \
        <(cat ~/data/5_pathways/gene_functions.txt)  \
        | cut -f2,3,4 \
        | sort -t$'\t' -g -k1,1 -S100M \
        | gzip \
        > ${E_DIR}/results_gene_functions_kegg.txt.gz

    ##
    ## functions_cnt<TAB>function
    ##
    zcat ${E_DIR}/results_gene_functions.txt.gz \
        | cut -f3 | sort -t$'\t' | uniq -c | sort -k1,1 -nr \
        > ${E_DIR}/results_gene_functions_uniq.txt

    ##
    ## functions_cnt<TAB>function
    ##
    zcat ${E_DIR}/results_gene_functions.txt.gz \
        | head -n1000 | cut -f1,3 | sort | uniq -c | sort -k2,2 -gs | sort -t$'\t' -k2,2 -su | sort -k1,1 -nr | sort -k2,2 -gs \
        > ${E_DIR}/results_gene_functions_best_uniq.txt

    ##
    ## functions_cnt<TAB>function
    ##
    zcat ${E_DIR}/results_gene_functions_kegg.txt.gz \
        | head -n1000 | cut -f1,3 | sort | uniq -c | sort -k2,2 -gs | sort -t$'\t' -k2,2 -su | sort -k1,1 -nr | sort -k2,2 -gs \
        > ${E_DIR}/results_gene_functions_kegg_best_uniq.txt


    ##
    ## probability<TAB><gene_kegg_id<TAB>path_id<TAB>pathname
    ##
    join -t$'\t' -1 2 \
        <(join -t$'\t' -2 2 \
            ${BASE_DATA_PATH}/5_pathways/kegg_pathway_map_full.txt \
            <(zcat ${E_DIR}/results_gene_functions.txt.gz | cut -f1,2 | sort -t$'\t' -k2,2 -S100M) \
        | sort -t$'\t' -k2,2 -S100M \
        ) \
        <(cat ${BASE_DATA_PATH}/5_pathways/pathways_full.txt | cut -f1,4) \
        | awk -F'\t' '{print $3"\t"$1"\t"$2"\t"$4}' \
        | sort -t$'\t' -g -k1,1 -S100M \
        | gzip \
        > ${E_DIR}/results_pathways.txt.gz     

    ## 
    ## probability<TAB>path_id<TAB>intersection<TAB>pathname
    ##    
#    join -t$'\t'  \
#        <(6_annotation/intersect_genes_avg \
#            -m ${BASE_DATA_PATH}/5_pathways/pathway_vectors_avg_${E_VALUE}_.txt \
#            -v $E_DIR/vectors.txt \
#            -g $GENOMES_NUM \
#            | cut -f2- | sort -t$'\t' -k1,1 -S100M) \
#        <(cat ~/data/5_pathways/pathways_full.txt | cut -f1,4) \
#        | awk -F'\t' '{print $2"\t"$1"\t"$3"\t"$4}' \
#        | sort -t$'\t' -g -k1,1 -S100M \
#        > $E_DIR/results_pathways_avg.txt
#        | sort -t$'\t' -k4,4 -S100M -u \
#        | sort -t$'\t' -g -k1,1 -S100M \

#    break    
#done

